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Abstract 

The free energy of a hydrogen-helium fluid mixture is evaluated for the 
temperatures and densities appropriate to the deep interior of a giant planet 
such as Jupiter. The electrons are assumed to be fully pressure-ionized 
and degenerate. In this regime, an appropriate first approximation to the 
ionic distribution functions can be found by assuming hard sphere inter- 
actions. Corrections to this approximation are incorporated by means of 
the perturbation theory of Anderson and Chandler. Approximations for the 
three-body interactions and the non-linear response of the electron gas to 
the ions are included! We predict that a hydrogen-helium mixture, containing 
10% by number of helium ions, separates into hydrogen-rich and helium-rich 
phases below about 8000°K, at the pressures relevant to Jupiter (4-40 
Megabars). We also predict that the alloy occupies less volume per ion 
than the separated phases. The equation of state and other thermodynamic 
derivatives are tabulated. The implications of these results are mentioned. 
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Introduction 

The mass and radius of a giant planet such as Jupiter can only be explained 

by assuming that the main constituent is hydrogen.* - This suggests that Jupiter 

may have roughly solar composition so that about one atom in ten is helium. More- 

2 

over, Jupiter emits about twice as much radiation as it receives from the sun. 

This indicates an internal heat source and is consistent with a temperature in the 

deep interior that exceeds the melting point of metallic hydrogen or helium.* - It 

has also been suggested that the helium may have only limited solubility in the 
3 4 

hydrogen. ’ Clearly, detailed models of the giant planets require an understanding 
of the thermodynamics and phase diagram of dense hydrogen-helium fluids. In this 
paper, the relevant properties of such a system are calculated. We assume that all 
the electrons are pressure-ionised, although our calculations have at least approxi- 
mate validity at the lower pressure relevant to Jupiter. 

The only previous relevant calculations are the incomplete Monte Carlo studies 
5 6 

by Hubbard. ’ The lengthy computational time of those calculations is avoided 

here by choosing an appropriate trial solution for the ionic distribution functions 

and then using a perturbation theory (the optimised random phase approximation of 

Anderson and Chandler^) to closely approximate the real ionic distribution. At the 

densities and temperatures of interest (1 ^ p ^ 10 g/cm , 10 < T ^ 5 x 10^ 0 K), 

the ion-ion interaction (in the presence of screening electrons) is characterized 

8 

by a repulsive core and a weak long-range part, so the appropriate trial solutions 
are the distribution functions for hard spheres. At higher densities or temperatures, 
a different approximation scheme (or the Monte Carlo method) is more appropriate. 

Our calculation contains features not present in Hubbard's calculations. We 
have used a more realistic dielectric function, and corrections have been made for 
the quantum mechanics of the ions, the three body interactions, the non-linear 
response of the electron gas to the ions, and finite temperature corrections to tte 
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electron gas. Some of the thermodynamic properties are significantly affected by 
these corrections. Nevertheless, our results are similar to those of Hubbard in 
most instances. We do make an important new prediction: the existence of a 

miscibility gap in the hydrogen-helium alloy. We also predict that there is a 
small but non-negligible volume non-additivity (i.e. the alloy occupies less volume 
per ion than the separated phases). The detailed astrophysical implications will 
be discussed elsewhere; they are only briefly mentioned in the present paper. 

In Section II, the free energy calculation is described. In Section III, the 
thermodynamic properties are discussed. In Section IV, the phase separation is 
described and in Section V, the validity of our calculation is assessed. 



II. The Free Energy 


Consider a binary system in which a fraction x of the nuclei have charge Z = 2, 

and a fraction 1 - x have charge Z - 1. The compensating electron gas has an 

4 3 3 -1 

average density Z*N = ( "3 rrr s a o ) where N is the ion number density, the effective 
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valence is Z* = 1 + x, a Q - 0.529 x 10 cm, and r g is the electron spacing parameter 

(r ~ 1 in Jupiter). The free energy is written as 
s 


F = F + F. + E + E + AF . . + F rt 
eg hs M BS mt Q 


( 1 ) 


We now discuss the meaning and evaluation of each term. 

The free energy of a uniform electron gas, F , can be expressed as the sum 
of a zero temperature contribution and a small finite temperature correction since 
we assume 

( 2 ) 


T < < 6 X 2 10 °K = Tp 


s 


where T and T^ are the actual and Fermi temperatures, respectively. The zero 
temperature contribution is 

E - (2.21/r 2 - 0.916/r - 0.115 + 0.0310m: )Z* Ryd/ion 

Gg S S S 


(3) 
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where an interpolation approximation has been used for the correlation energy. 

(Any small inaccuracy in the correlation energy has negligible effect on the final 

results). The only finite temperature correction that we have retained is the 

10 


kinetic energy term 


2 2/3 
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eg 


ion 


( 4 ) 
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where k fi T is in Rydbergs (k g = 6.34 x 10 Ryd/°K). The other corrections can be 
shown to be negligible at the temperatures of interest. 

As a first approximation, we assume that the ions, with their comoving screening 
clouds of electrons, interact like hard spheres. Consequently, we have, in addition 
to electrostatic energy contributions, the free energy, F hg , of an equivalent neutral 
hard sphere liquid. This has been approximately calculated by Mansoori et all -2 : 
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where y^ = x(l - x) (1 - a) (1 + a)/d 

y 2 = x(l - X)a(l - a) 2 [(l - x)a 2 + x]/d 2 
y 3 = [(1 - x)a 2 + x] 3 /d 2 

3 

d - (1 - x)a + x 

- volume occupied by hard spheres 
^ ~ total volume 

diameter of hydrogen hard spheres 
a “ diameter of helium hard spheres 

and are the hydrogen and helium ionic masses respectively. is the free 

energy of an ideal gas mixture. (A hard sphere mixture deviates from an ideal gas 

to the- extent that tj is non-zero) . For r) ^ 0.45, a classical liquid is expected 

13 

to solidify. The predictions of equation (5) are almost indistinguishable from 

hard sphere Monte Carlo calculations, in contrast to the alternative formula of 

14 

Lebowitz. The parameters r\ and a are not known in advance and are to be determined 
variationally . 

We next evaluate the Madelung energy E^, which is the electrostatic energy of 
point ions immersed in a uniform electron gas. 
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where the integrals extend over the entire volume Q and (...) denotes an ensemble 
average. In ( 6 ), r^ refers to a nucleus of charge Z ± at position r n . 


It can then 
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be shown that 
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( 7 ) 


(Z^ = 1, Zj * 2 for H-He), where the partial structure factors S^(k) are defined as 

,, v 1 / V ik* (r* - - rb \ 

s ( k ) = — < ) e ~ v n ^tn ) 
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and fL,TL are the partial number densities of the two ionic species. We approximate 
by the Percus-Yevick hard sphere result found by Ashcroft and Langreth.^ Ross 
and Seale* have shown that in this approximation, the summation in (7) can be 
evaluated exactly. We give the rather complicated result here because of an error 
in their paper: 




-JH 
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ion 
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( 9 ) 


2 10 2 2 
G n = [Tjj - 2^ + 10 - ^T| 1 T 1 2 + 10ri 2 + — Cti 1 T| 2 - n 2 ) + (ri 2 - 2 T 1 2 - 411^) /a ] 

G 12 = 15(T1 2 ■ V ( a " L >2a/(l + a > - 30 < 1 + 'n) 

+ (2a '1 + a) 2 [rt 2 - 2^ + 10 - 4-rj^ + 20ri 2 + 10 (r^ + r\ + 2) /a 
+ (ri 2 - 2 ti 2 + 10 - 4ti 1 ti 2 + 20 rip/a 2 ] 

3 

= (1 - x)a r)/d, r\ 2 = xri/d 

and G 22 is obtained from G^ by replacing - by a and by r^. Corrections to the 

hard sphere approximation are contained in AF^ , discussed below. 

We next consider the "band structure" energy E , 

BS 


resulting from the non- 
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uniformity of the electron gas. In general, it is possible to expand 

the Fourier transform of the electron density change induced by the ions of species 

17-19 

i, in powers of the electron-ion interaction 

= 1 oJS’os) 

n-1 
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We have retained terms up to third order in the electron-ion interaction, so 

=» r.(2)j „(3) 


BS 


'+ E 


where 


E<n) ’ din Z < Z <-S> > 


(12) 


MO i, j-1,2 

Note that in equations (7) and (12), the zeroth Fourier components are omitted 
from the summation since they are exactly compensated by the divergent electron- 
electron interaction energy. 

The lowest order term in p* ,(k) is the linear response result^’ ^ 


- iff ( efe) - «9 


k 2 ( 1 


(13) 


where e(k) is the static dielectric function. In our calculations, the Hubbard 

20 

approximation has been used. This is a more realistic form for the linear response 

than the Lindhard expression used in the Monte Carlo studies'*' Corrections for 

21 

the dynamic response of the electron gas appear to be negligible and have not 


been included. 



- 8 - 


that 


From the definitions of the partial structure factors, it can then be shown 
22 


E<2) = ^ J { z i^ " x ) s u 00 + 2 Z x Z 2 x^( 1 - x) % S 12 (k) + Z 2 xS 22 (k)} 

X leTk) " R y d/iori 


-1 


where k is in units of a . This contribution is evaluated numerically using the 

o 

hard sphere structure factors. 

The theory for is not given in detail here, since it is an obvious 

17-19 

generalization of the results for crystalline metals . (For example, the only 

change to equation 4.8 of Lloyd and Sholl^ is the insertion of ensemble averages.) 

19 

The appropriate generalization of equation (90) of Hammerberg and Ashcroft is 

1/3 


,<3) _ 
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where T 




(16) 


m,n,p 


W 1 (§) = — , ? = 77 ” , k is the Fermi wavevector, and ’ is the three- 

H(5) V 

vertex function defined by Hammerberg and Ashcroft. Since m,n,p do not have to be 

(3) 

different in the summation for T„^> it is clear that E contains not only 

three-body interactions, but also corrections to two body interactions and (structure 

independent) corrections to the screening of single ions. The evaluation of T^.^ 

can only be approximate since there is no accurate theory for three-body correlation 

functions. We have chosen a convolution approximation (see Appendix I). 

We have not evaluated the fourth order contribution to the band structure 

energy, but some semi-quantitative assertions can be made. First, it is clear from 

19 

the work of Hammerberg and Ashcroft that there are additional complications at 


fourth order that can only be encompassed by the use of finite temperature perturba- 
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tion theory. Second, their formal results can be easily extended to liquids and 

liquid mixtures. In particular, the terms which they ascribe to the non-sphericity 

of the Fermi surface in the solid, are non-zero (and comparable in magnitude) in 

the liquid phase. It seems that the substantial cancellation of fourth order terms 

that they found for crystalline metallic hydrogen, persists in the liquid. The 

cancellation is also substantial for helium, since the dominant fourth order con- 

(4) 

tributions have similar Z dependence. Thus, it is hoped that our omission of E 
is not a serious deficiency. 

In the above calculations, we have used an effective ion- ion interaction, the 
accuracy of which is limited only by uncertainties in the dielectric function of 
the electron gas. We have not, however, evaluated the ionic configuration appropriate 
to that interaction. This is corrected by the optimised random phase approximation 
method of Andersen and Chandler(see Appendix II). To second order in the electron- 
ion interaction, 

AF. nf . « -^ 3 - J W F det l* + ^ /k R T l 1 " mi & ^ /k R T I d ls d 7 ) 

in 16 tt N J 


where "det" means determinant, "Tr" means trace, I is the unit 2x2 matrix; and 
sT, $ are the 2x2 matrices, the elements of which are 

r ij - V*> 


. . = (x.x.)^cp. . (k) 
lj 1 J ~ 


(18) 


where x^, are the number fractions of ion species j and j; and cp„(k) is the 
Fourier transform of the optimised potential (r) , given by 


U t .(r) = - V^r) 


= vef f (r) 
ij ' ' 


r < R. . 
ij 


r > R. . 

“ 13 


(19) 


eff . 


where R. . is the minimum hard sphere separation and v . (r) is the effective ion-ion 

J 

2 2 

interaction, the Fourier transform of which is 4rre /k e(k). The functions C„(r) 


are chosen variational ly 
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S ^ F int 

ficTTc?) ■ r < R ij 


A quadratic function of r was found to be adequate for approximating each C.^(r). 

Finally, we have included the lowest order high temperature quantum correction 


F = JOL 

Q 24k B T 


fg. .(r) V 2 v Cff (r)dr 
M, . J 11 n w ^ 


i. J-1.2 


ij i J 

where g„(r) is the pair correlation function for ion species i and j, and FL is 

23 

the mass of ion species i. This result was first derived by Wigner and is the 
effect of including the uncertainty principle for the ions, to lowest order in 
h /Ma o kgT . Unlike F^, the next term in the Wigner expansion depends on whether 

.O 

the ions are bosons or fermions, but it can be shown to be negligible (~ 10 Ryd) 
for our present purposes. 

-3/2 

Notice that F does not scale as r the result that one might expect for 

x S 

24 

zero -point motion. The Wigner expansion is rapidly convergent provided 


Mo V 


< < i 


where a is roughly the range of the strongly repulsive part of the effective inter- 
action. Detailed calculations indicate r T ^ 1500°K is a sufficient condition for 

s 

Fq to represent accurately the quantum correction. 



III. Thermodynamic Properties 


The free energy given by equation (1) was evaluated as a function of density, 

temperature and relative concentration. The procedure is to minimize F^ g + E^ 

( 2 ~) ( 3 \ 

+ E v ' + J + Fq as a function of T| and a. (The best choice for a is found to 

be insensitive to all other parameters and varied between 0.73 and 0.77 only). The 

remaining contributions to F are then added. Since AE . effectively corrects the 

xnt 

hard sphere approximation, to second order in the electron-ion interaction, the 

( 3 ) ( 3 ) 

minimization procedure would be unnecessary, were it not for E . If E is 
excluded, then the total F is indeed very weakly dependent on T| and a and this 
encourages confidence in our calculation. The minimization procedure is justified 
by a result of thermodynamic perturbation theory which states that the exact free 
energy is bounded above by the free energy calculated using the hard sphere model. 

Rather than tabulate F, we have tabulated various derivatives that are par- 
ticularly useful in constructing planetary models(see Tables I-III) . Note that 
the heat capacity per nucleus, at constant volume, can be determined from Tables 
I and II since 


25 


c 
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where 


= — f- ) 
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S ,x 


The equation of state P(V,T,x) has not been tabulated but has been fitted to a 

polynomial in x and r for T = 6000°K. (Pressures at other temperatures can be 

s 

found by using Table 1). 


p = [1 + a(x)r g + b(x)r^ 4- c(x)r^]Mbars 

r 

s 


(24) 
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a(x) = -0.654 - 0.200x + x(l - x) (-0.182 + 0.370x - 0.288x ) 

b(x) = 0.085 - 0. 054x + x(l - x) (-0.086 - 0.530x + 0.573x 2 ) 

c (x) = -0.008 + 0.028x + x(l - x) ( 0.077 + 0.254x - 0.321x 2 ) 

This interpolation formula is accurate to 0.37, for 0.6 ^ r ^ 1«3» For x = 0 

s 

(metallic hydrogen) our results for P are in good agreement with previous workers. 

5 

For example, at r g = 1.0 we obtained P= 21.8 Mbars whereas Hubbard and Slattery 

obtained 22.0 Mbars. This agreement is not surprising, since the pressure is 

determined mainly by E and E^, terms that are common to both treatments. (For 

example, contributes only 0.35/r 2 Mbar to P for hydrogen). We believe that 

s 

our equation of state is the best available for liquid metallic hydrogen and is 
probably accurate to better than 1%. The accuracy for x 4 0 is more difficult 
to assess, since the perturbation expansion is much less valid for helium (as we 
discuss in the final section.) 

The parameter c (Table I) would be 1,0 for an ideal gas and 1.5 for a high 
temperature Debye solid. The actual behavior is more complicated than either of 
these limiting cases. For example, c is reduced by the quantum corrections at low 
T but increased by the free energy of the electron gas at high T. The Monte Carlo 
results'* are too incomplete for a detailed comparison. To give a sample comparison: 
at x = 0, r g = 1.0; they obtained 1.49 at 4200°K, 1.42 at 6300°K, 1.32 at 10500°K 
and 1.16 at 31500°K. Our results for the same temperatures are 1.04, 1.15, 1.18 
and 1.21 respectively. 

The parameter y(Table II) would be 2/3 for an ideal monatomic gas and 0.5 for 

a high temperature Debye solid(neglecting screening). The actual value usually lies 

between these limiting cases. For comparison, the Monte Carlo result for x = 0, 

3 4 

r = 1.0 was y ~ 0.64, for 4 x 10 ^ T ^ 3 x 10 °K. 

s 

The specific entropy (Table III) differs surprisingly little from the Monte 

3 26 

Carlo results. For example, at x = 0.143, p = 5 g/cm Hubbard finds T = 23'400°K 
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for S = 8.4 k^/nucleus. Our results predict T = 20 , 500°K. (For a fully adiabatic, 
homogeneous Jupiter, this would be roughly the central temperature of the planet.) 

Our calculations enable us to assess accurately the "volume additivity" 
approximation that has often been made. Let fi(x,P) be the volume per ion. We 
can always write 

[1 + 6(x,P) ]Q(x,P) = xO(l,P) + (1 - x)fi(0, P) (25) 

where 6(x,P) is a measure of the deviation from volume additivity in the alloy. 

Figure 1 shows that 6 is significantly non-zero, in contrast to Thomas-Fermi 

27 

theory where it is natural to assume 6=0. This non-additivity is not attributable 
to any particular term in the free energy. It is comparable to (but usually 

2 

somewhat larger than) the non-additivity observed in liquid alloys in the laboratory. 

This result indicates a small modification to models for giant planets. For 
example, a model constructed using the exact equation of state and x = 0.1, would 
require x =** 0.12 if volume additivity were assumed. Correct allowance for non- 
additivity slightly reduces the amount of helium required in giant planet models. 



IV. Phase Separation 


There are two ways of testing for incomplete miscibility in a fluid mixture 

calculation. One way is to look for divergent behavior in the long wavelength 

limit of a partial structure factor, corresponding to the onset of macroscopic 

29 

concentration fluctuations. This method has been applied by Stroud to the alkali 
metals with some success, but his mean field approach would predict complete misci- 
bility in the H-He system. The more exact calculation described below indicates 
that this must be a failure of the mean field approximation. (Our inability to 
find divergent behavior indicates a lack of self-consistency in our calculation. 

It is hoped that this inadequacy is serious only near the critical point.) 

The second test is to evaluate the Gibbs energy G(P,T,x) and then determine 

whether & 2 G/dx 2 < 0 in any region of (P,T,x)-space. Such regions are unstable 

30 2 2 

towards phase separation. This method has been used to predict, with considerable 
success, the miscibilities and phase separation curves for several alkali metal 
mixtures . 

Figure 2 shows the Gibbs energy of mixing, defined as 

AG(P T, x) = G(P,T,x) - xG(P ,T , 1) - (1 - x)G(P,T,0) (26) 

At low temperatures, the unstable region is easily discemable. Near the critical 
temperature, a careful common tangent construction must be made (see the curve at 
9000°K, for example). Since AG is much smaller than G, it is clear that even small 
errors can dramatically affect an estimate of T , the critical temperature. Never- 
theless, there is little doubt that the unstable region exists for T ^ 7000°K. 

In Figure 3, we show the phase separation curves constructed from several plots 
like Figure 2. (Not shown are the results at P = 200 Mbars, for which T £ 10'000 a K) 
The regions near T c are interpolations and may be inaccurate. Away from T =“ T c> the 
curves are likely to be accurate to about + 20%. The results may also be incorrect 
at low temperatures (T < 2000 e K) where solid phases may exist. We found no evidence 
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in the calculations for the more complicated phase diagrams that are permitted by 

31 

the Gibbs phase rule. 

The results indicate that a H-He mixture of solar composition (x 0.1) 

separates into hydrogen-rich and helium-rich phases at temperatures less than about 

8000°K and pressures in the range 4-200 Mbars. It is not surprising that this 

phase transition did not clearly manifest itself in Hubbard's Monte Carlo calculations, 

32 

because of the small number of particles he used. The immiscibility is not 

attributable to any particular contribution to F. 

The effect of separation in a planet such as Jupiter is to retard the cooling 
t\ 

rate and evolution. This will be discussed elsewhere. 



V. Discussion 


We have assumed throughout that it is valid to consider the helium as fully 
pressure-ionized, even at a few megabars pressure. There may be serious objections 
to this. 

It should be emphasized, however , that the validity of our approach has nothing 

to do with the pressure at which pure solid helium becomes metallic. We have 

evaluated the band structure of face-centered cubic helium using plane waves as a 

basis set for the electronic wavefunct ions . We obtain a transition pressure of 

33 

70 Mbars, similar to the result of Trubitsyn „ However, a calculation of the 
band gaps to third order in the electron-ion interaction is accurate even at 10 Mbar. 
This is a more relevant criterion, since our calculation of F relies on the con- 
vergence of a plane wave expansion and not on the existence of metallic conduction 
in the helium fluid. We have also used the methods outlined in this paper to cal- 
culate the free energy of molecular hydrogen at r ^ 1.7 to an accuracy comparable 

s 

to that achieved using semi-empirical pair potentials. We mention this to 

emphasize the power of the perturbation technique used. 

3 

Contrary to what has been stated in the literature, there is no good reason 
for supposing that helium becomes less soluble at pressures lower than those 
considered in this paper. Indeed, the screened interaction between a neutral 

34 

helium atom and a proton has large attractive region, which suggests miscibility. 

In contrast, the solubility of helium in alkali metals at near-zero pressure is 

35 36 

very low because of the repulsive electron-helium pseudopotential . As the 

pressure increases and the wavelength of a conduction electron becomes comparable 

to the size of a helium atom, the pseudopotential becomes less repulsive and the 

solubility increases. It is not correct to compare the H-He system with any large 

r system (such as Na-He) that is accessible in the laboratory, 
s 

In this paper, we have shown how a judicious mixture of perturbation techniques 
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enables us to evaluate the thermodynamics of a non-ideal system that was previously 
thought to require Monte Carlo techniques. It is likely that other, similar systems 
will yield to a comparable analysis. One candidate is the ^-He mixture that is 
present at lower pressures (P ^ 2 Mbars) in the giant planets. 
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Appendix I 


Evaluation of T 


lJk ( Sl-3Z> 


The calculation of E 


(3) 


requires an approximation for 




< -2 

NO 


l ' 

m, n,p 


+ ( 5i ' + 


(Al-1) 


For simplicity, we consider a pure liquid so that the subscripts i,j,k can be 

37 38 

omitted. There are numerous papers in the theory of liquid metals ’ in which 
T(<Ji>^ 2 ) i- s approximated by T c (q^,q 9 ), where 


T c <*!»&> = s <“3l )S( S.l - l2 )S( & 


(Al-2) 


and S(q) is the usual liquid structure factor. This is often called the "geometric 

37 

approximation". As discussed by Ballentine and Heine , it treats clusters of 

three atoms approximately, but is otherwise exact. What has apparently not been 

pointed out before is that (Al-2) is identical to the well-known convo lution 

39 

approximation in real space. This approximation states that 


g(3) (£l^2^3 ) “ ' (8 <2) <^l’l2 ) “ 1 Ke (2> <^2’-3 ) ‘ 

x (S (2) <£ r r 3 ) - 1) + N J(g (2) (r 1 ,r 4 ) - 1) (g (2) ( rj - 1) 


.(2), 


x <8 ( £ 3 ^ 4 } ' 1)d U 


(Al-3) 


(2) (3) 

where g and g are the two-body and three-body correlation functions respectively, 
(3) 

and is the Kirkwood superposition approximation 


®KS ^r~2’^ ~ g ^£ 2 ’ *3^ 


(Al-4) 


To prove the equivalence of (Al-2) and (Al-3) we note that by definition 


/ 1 V i(-q-*r + (q. - q_)*r + q„* r )\ 

'NO Zj 6 Mn -W- ^-2 ->-n -J2 V> 


m,n,p 
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= J dr ds dt <4l - S2'tW 3) (j^s,t) (Al-5) 

^ * 

where ) means m, n and p are all unequal in the summation. The proof is then 
straightforward. 

39 

As discussed by Feenberg , the convolution approximation is an exact solution 

(2) (3) * - 

of the hierarchy equation that links g v and , whereas the superposition 

approximation is not. Moreover, the convolution approximation is a natural conse- 

40 , s 

quence of diagrammatic analyses, such as that made by Abe. However, it does not 

f 3) 

necessarily satisfy the physical requirement, g v > 0, whereas the superposition 
approximation does. This could lead to serious errors for strongly repulsive 
potentials. 

Nevertheless, we have used the convolution approximation since the super- 
position approximation is very cumbersome in Fourier space: 


T Ks ( ar*s> = ' (s( -Si> - D^i ' S2> - ~ 

dk (S(k) - 1) (S (k + %) - l)(S(k + ^) - 1) 


+ -V-J 

(2tt) 3 N J 


(Al-6) 


(3) 

We have made one test evaluation of E v using instead of T^, for T| = 0.3 

and r g = 1.0 in pure hydrogen. The results agreed to within 10%, although T rS 

and T often differed by more than 10%. This is also expected to be comparable to 
c 

the error that is incurred when either approximation is used, rather than the "true" 
T(qi,q 2 ) that would be determined by Monte Carlo or molecular dynamics techniques. 

For larger T| (closer packing), the error may be larger, since it is observed in 

(3) (3) 41 42 

machine calculations that g^' deviates more from the true g as T) increases. ’ 

The generalisation of (Al-2) to mixtures is straightforward (but not trivial). 

The result is 

T ijk ( ai-S2 ) ■ (S ij ( ‘Jl ) ' V<VSl ‘ ‘ V <S ki < 32 ) ‘ 6 ki } 

4,0 '-'a 


+ <Vi> (S 13 (-5i) - V ( VSi - v 


V + (x jV < s tj ( -s.i> - 6 ij 
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X (S ki ( V “ 6 ki> + (x i X j ) (S jk% " S2> " 6 jk )(S ki ( -S 2 ) - 6 ki> 

+ 6ij < x I * b > %S lk (52> + 6 tk<Vj>Slc<ai + 6 jk (X i X j ) ^ij ( 'ai ) 


- 2 *l 6 ij 6 jk 


where 6.. is the Kronecker delta, 
i J 


(Al-7) 



Appendix II 


The Optimised Random Phase Approximation (ORPA) 


The object of a liquid perturbation theory is to approximate the true pair 
interaction by a very simple interaction for which the corresponding ionic con- 
figuration is well known. The free energy is then expanded about the free energy 
of the well understood reference system, in powers of the difference between the 

7 43 

actual and reference interactions. The ORPA was devised by Andersen and Chandler ’ 

so that this perturbation expansion would rapidly converge. In briefly discussing 

our application of this method, we restrict ourselves to a pure system. 

eff 

The first step is to approximate the real interaction v (r) by a trial 
interaction v T (r) given by 

v T (r) = » r < R 

= v e ^(r) r > R (A2-1) 

where R is the variationally determined hard sphere diameter. The error in F 

that is Incurred by thi'is replacement is 

Nk B T “ eff 

= Hr J C x C£> ex Pt- v (£> /k B T J d £ (A2-2) 

r<R 

where C,j,(r) is the direct correlation function, evaluated by equation (A2-5) below, 
eff 

Since v (r) > > k^T when r < R, this error is found to be negligible (it corresponds 
to neglecting infrequent high-energy collisions between ions). 

The trial interaction is then decomposed into a reference part and a pertur- 
bation part. 

v T (r) = V Q (r) + u(r) 

v (r) = » r < R 

= 0 r > R 

u(r) = -k fi Tc(r) r < R 
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An exact calculation of F for this interaction roust be independent of c(r), since 
it is defined in the physically inaccessible region. However, the result of an 
approximate calculation does depend on c(r), Andersen and Chandler showed that 
AF int , the change in the free energy from a hard sphere system, can be accurately 
approximated by the lowest order term 


AF 


( 1 ) _ 

int 


V 

16rr\ 


J> [1 + Hgaafa ] . sma} («.*, 


where cp(k) 


-J 


ilc* IT 

e ^ ~ u(r)dr, provided 


6AF 


( 1 ) 

int 


6c(r) 


= 0, r < R 


(A2-5) 


They also showed that c(r) C_,(r) - C, (r) for r < R, where C, is the hard 

i ns hs 

sphere direct correlation function. 

In our calculations j AF iTlt | ^ kgT, but is large enough to have a significant 
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effect on the entropy, as the Monte Carlo calculations by Ross would have pre- 
dicted. 

In Figure 4, a comparison is made of the pair correlations for ORPA, hard 
spheres and Monte Carlo 5 . The similarity of our ORPA result and the machine 
simulation is very satisfactory. 
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TABLE I 


c 




V 


r 

s 

r T 

s 




X 




0.0 

0.1 

0.2 

0.4 

0.6 

0,8 

1.0 

r.2 

4000 

1.07 

1.12 

1.18 

1.40 

1.60 

1.46 

1.30 


7000 

1.24 

1.30 

1.41 

1.67 

1.86 

1,78 

1.70 


12000 

1.22 

1.25 

1.30 

1.39 

1.52 

1.69 

1.83 


20000 

1.23 

1.25 

1.27 

1.32 

1.41 

1.47 

1.55 


30000 

1.25 

1.30 

1.31 

1.35 

1.44 

1.47 

1.55 

0.9 

4000 

1.01 

1.04 

1.05 

1.02 

0.83 

0.77 

0,69 


7000 

1.12 

1.16 

1.19 

1.32 

1.43 

1.51 

1.48 


12000 

1.17 

1.23 

1.28 

1.46 

1.56 

1.59 

1.65 


20000 

1.17 

1.19 

1.24 

1.32 

1.40 

1.49 

1.60 


30000 

1.19 

1.21 

1.25 

1.33 

1.40 

1.45 

1.59 

0.7 

4000 i 

0.65 

0.60 

0.57 

0.59 

i 

0.61 

0.64 

0.65 


7000 

1.03 

1.05 

1.10 

1.12 

1.12 

1.09 

1.08 


12000 

1.14 

1.18 

1.25 

1.37 

1.48 

1.54 

1.56 


20000 

1.18 

1.24 

1.30 

1.36 

1.41 

1.47 

1.50 


30000 

1.20 

1.25 

1.30 

1.35 

1.40 

1.44 

1.47 
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TABLE II 


Y = (hS?) s 


ir 

s 

r T 
s 

X 


0.0 

0.1 

0.2 

0.4 

0.6 

0.8 

1.0 

1.2 

4000 

0.61 

0.62 

0.62 

0.65 

0.66 

0.61 

0.55 


7000 

0.61 

0.62 

0.66 

0.67 

0.68 

0.61 

0.56 


12000 

0.61 

0.61 

0.65 

0.65 

0.66 

0.65 

0.65 


20000 

0.63 

0.64 

0.64 

0.65 

0.66 

0.65 

0.65 


30000 

0.64 

0.64 

0.65 

0.65 

0.66 

0.65 

0.65 

0.9 

4000 

0.58 

0.58 

0.59 

0.58 

0.55 

0.53 

0.53 


7000 

0.60 

0.61 

0.60 

0.59 

0.60 

0.61 

0.61 


12000 

0.63 

0.63 

0.62 

0.62 

0.62 

0.63 

0.62 


20000 

0.64 

0.63 

0.63 

0.63 

0.64 

0.63 

0.63 


30000 

0.64 

0.63 

0.63 

0.63 

0.64 

0.63 

0.63 

0.7 

4000 

0.56 ' 

0.56 

0.55 

0.55 

0.55 

0.53 

0.53 


7000 

0.59 

0.60 

0.59 

0.59 

0.59 

0.56 

0.56 


12000 

0.63 

0.62 

0.62 

0.62 

0.62 

0.62 

0.63 


20000 

0.63 

0.64 

0.63 

0.63 

0.63 

0.63 

0.64 


30000 

0.63 

0.64 

0.64 

0.63 

0.64 

0.63 

0.63 


R 
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TABLE III Entropy S (in k /nucleus) at r =0.9. 

B s 

(For extrapolation to other densities, 
use Table II.) 


T(°K) 

X 

0.0 

0.1 

0.2 

0.4 

0.6 

0.8 

1.0 

3500 

4.3 

4.8 

5.15 

5.6 

5.9 

6.0 

5.75 

6000 

5 ♦ 25 

5.75 

6.1 

6.5 

6.75 

6.85 

6.55 

9000 

6.0 

6.55 

6.95 

7.45 

7.75 

7.85 

7.55 

12500 

6.6 

7.2 

7.6 

8.2 

8.6 

8.75 

8.6 

17500 

7.25 

7.85 

8.3 

8.9 

9.4 

9,6 

9.5 

25000 

7.95 

8.55 

9.0 

9.65 

10.15 

10.45 

10.45 

35000 

8.55 

9.2 

9.65 

10.35 

10.85 

11.25 

11.25 
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Figure Captions 


4 

Figure 1 Extent of Volume Non-Additivity (equation (25)) at T = 10 °K and for 
various pressures. 

Figure 2 Gibbs energy of mixing AG, as a function of T and x at P = 8Mbars. 
Figure 3 Phase separation curves for various pressures. The phase-excluded 
region (miscibility gap) is below the curve in each case. 

Figure 4 Pair correlation functions for Z = 1, r g = 1.0 and T = 4200°K. 

=— — — optimised hard spheres 

ORPA 

5 

Monte Carlo (Hubbard ) . 


29 












